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CLASSICAL EIGHTH- AND LOWER-ORDER 
RUNGE-KUTTA-NYSTROM FORMULAS WITH STEPS IZE CONTROL 
FOR SPECIAL SECOND-ORDER DIFFERENTIAL EQUATIONS 


INTRODUCTION 


1. In earlier reports [ 1J , L 2] , [3] the author derived and listed Runge-Kutta 
formulas with stepsize control for first-order differential equations. 
Although second-order differential equations can be converted to first- 
order differential equations by introducing the first derivatives as new 
variables, this amounts to an increased computational effort, since we 
then have to deal with twice as many differential equations as in the origi- 
nal second-order problem. Therefore, the direct Runge-Kutta integration 
of second-order differential equations without conversion to first-order 
equations might be preferable. Direct Runge-Kutta formulas for second- 
order differential equations were first published by E. J. NYSTROM [ 4] . 

We therefore refer in this report to such direct Runge-Kutta formulas for 
second-order differential equations as Runge-Kutta-Nystr5m (RKN) 
formulas. 

2. In this report we restrict ourselves to a special class of second-order 
(vector) differential equations, 

x f(t,x) , (1) 

which do not contain the first derivative x on the right-hand side. Such 
special second-order differential equations are frequently encountered in 
mechanics and physics. 

The derivation of the equations of condition for the Runge-Kutta -Nystrom 
coefficients is much simpler and easier for the special equations (1) than 
for the general equations which would also contain the first derivatives 
x on the right-hand side. 

3. Similar to the Runge-Kutta formulas of our earlier reports [ 1] , [2] , [ 3] , 
the Runge-Kutta-Ny strom formulas of this report include an automatic 
stepsize control based on a complete coverage of the leading term of the 
local truncation error in x. This coverage is achieved by one additional 
evaluation of the differential equations. Each of our Runge-Kutta-Ny strom 
formulas represents in fact a pair of integration formulas for x which differ 



from one another by the one additional evaluation of the differential equa- 
tions. The orders of these two formulas differ by 1. Therefore, the 
difference of the formulas represents an approximation for the leading 
term of the truncation error in x for the lower-order formula. By 
requiring that this difference remain between preset limits, an auto- 
matic stepsize control for the lower-order formula can be established. 


SECTION I. EIGHTH-ORDER FORMULA RKN 8(9) 


For the sake of brevity in the derivation of the formulas, let us consider 
in the following a system of two second-order differential equations: 


x - f (t,x,y) 
y = g(t,x,y) 


( 2 ) 


Our results, however, will hold, in a quite obvious way, for a system of 
any number of second-order differential equations. We introduce 


f 0 = f(t 0 , x 0 , y 0 ) 

/ K-l 

f « ‘ + ■I'V.h'l' • Z 

\ A-0 

K - 1 

y 0 + y 0 « t ,h + h 2 - V y g 
h A^O KA X 

(k - 1,2, . . . , 11) 

and corresponding expressions for the second equation (2). 


(3) 


For any k the sum of the coefficients y (A = 0, 1, .... k - 1) is re- 

K A 

lated to a by: 
k J 


K-l 

z 


A=0 


y 


kA 


J. 

2 


a 


(4) 


as can easily be seen by a Taylor expansion of the x- or y- argument of 
f in (3). 



For an eighth-order Runge-Kutta-Nystrdm formula with stepsize con- 


trol we then require 

10 

x = x 0 +x 0 h+h 2 • Yj c K f K +0(h9) 

K=0 

11 

^ = x 0 +x 0 h+h 2 • Y ^ f + ° (hl0) 
K=0 KK 

10 

x = x 0 + h • V c f + 0(h 9 ) 

u LJ „ K K 


15) 


with 

c = c for k - 0, 1, 2, ... , 9 

K K 


A 

c 10~ 


0 


A 

C H= c 10 


(6) 


and corresponding formulas for the solution of the second differential 
equation (2). In (3) and (5) the quantities to, x 0 , y 0 , x 0 » Yo are the 
initial values for the integration step under consideration, while h 
stands for the integration stepsize. 

As in our earlier report [2] , we require that the last evaluation of the 
differential equation can be taken over as first evaluation for the next 
step, thereby reducing the number of evaluations per step by one. By 
this requirement, the coefficients an d a ^ are determined: 


Ylio ~ c o> Ym — c i> Yii 2 — c 2 > • • • » Ymo Cio, 1 • 17) 

Our problem then consists in finding the Runge-Kutta-Nystrom coef- 
ficients a , y , c , c so that the right-hand sides of (5) are really 

K K A K K 

eighth- or ninth- order approximations of x or x. By expanding in Taylor 
series the solution of (2) as well as the right-hand sides of (5) and 


3 



equating the corresponding terms in both series, equations of condition 
for these coefficients can be obtained. 


We first expand the solution x, y of (2) in a (truncated) Taylor series: 


10 x 


</') 


y 


o 


x = x ° + ■-> 
l>-i 

10 v 


y = yo + 


V 

V=1 


- h‘ 

iv) 

[) 

v\ 


( 8 ) 


using equations (2) and their derivatives for the computation of the total 
derivatives x 0 , y 0 in (8). 


For the computation of tlie derivatives of (2) we introduce the differential 
operator: 


8 .8 .8 

D = — + x — + v — 

at ax y ay 


(9) 


Obviously, the following rules hold for this operator: 

D (<P + ip) = D(< p) D(i p) 

D (cp • ip) - <pD(tp) + ipD ((p) ( (10) 


D[D n (</>)] = D n+1 (<^) + n 


D n 1 ( (p ) f + D n 1 ( ) g 
x y 


The operators D 2 (</>) , etc. , are defined as symbolic powers of D: 


D 


\ip) = (</? ^ (p • x + ip • y) 2 = <p + 2<p, • x + 2<p. • y 

\ t x y / tt tx ty J 


+ ip • x 2 + 2 <p • xy + W " 2 

XX vwr 


etc. 


xy 


yy 


Observing the above rules, the total derivatives of x are obtainable 
from (2). The resulting lengthy expressions for these derivatives are 
somewhat shortened by the introduction of the following abbreviations: 



( 11 ) 


{</>(f ) !/(f)} = <fi(f ) ^ ) >Mg> 

x x y 


<<p(f ) M f ) M f )> = M f ) M f ) + < C (f xy ) l M f ) ^2 (g) 

+ <Mg) + ^ ^ fyy ) ^l^g) 


( 12 ) 


Extending (11) to the case that ^(f ) is a product of two operators, 
iMf) = Mf)}, we define, 

{ V (n {*,<f x > ♦,(«}} = fiywyiitf] + *,« y > ♦• <8)l 

+ <p( f ) [i^l(g ) M f ) + Mg ) Mg)) 

y x y 


and similarly for the case that ip( f) in (11) is a product of more than two 
operators. 

In an obvious way, ( 12) can also be extended to the case of third- and 
higher-order partial derivatives f xxx > f xxxx ’ etc ‘ In the case of a third ‘ 

order partial derivative we define: 

o(f ) iMf) M f > M*)> 

XXX 

-<p ( f ) i^i(f) Mf) M^ 

XXX 

+ w(i ) li/ii(f) M f > Mg) + M f ) Mg) M f ) + Mg) M f ) M f >) 

xxy 

+ w(f ) IMf) Mg) Mg) + Mg) M f > Mg) + Mg) Mg) M*)l 

xyy 

+ <^ (f yyy ) Mg) Mg) Mg) 

and similarly for higher-order partial derivatives. 

Using the abbreviations ( 11) , ( 12) , we obtain from (2) for the total 
derivatives of x (always taken for t = t 0 ): 
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X II -f A 

X III =* D(f) 

x 1V - I) 2 ( f ) + {f f} 

X 

x v •-= r?(i) + 3 {n(f )f} + {f d (f)} 

X X 

X VI ■= iy*(f) + 6 ( n 2 (f )f} + 3 <f f*> + 4 { 0(f ) D(f )} 

X XX X 


f (f)} +{f {f f}} 


X X 


VII 


D 5 (f) + 10 { If^f )f} + 15 < D (f )f 2 > + 10 { D 2 (f ) D(f )} 


XX 


VIII 


+ 10 < f f D ( f ) > + 5 { D(f ) D 2 (f )} + 5 { D(f ) (f f}} 

XX X XX 

+ {f O 3 (f )} + 3 { f { D(f )f}} + {f {f D(±)}> 

X XX XX 

n 6 (f) + 15 {D 4 (f ) f} + 45 < (f ) f 2 >+ 20 {d 3 <f ) D(f )} 

X XX X 


+ 15 <f f 3 > + 60 < D (f )f D(f)> + 15 { D 2 (f ) rf(f)> 

XXX XX x 

+ 15 { O 2 (f ) {f f>} + 10 <f [13(f)] 2 > + 15 <f i 13 2 ( f )> 

XX XX XX 

+ 15 <f f{f f}> + 6 { D(f ) D?(f)} + 18 ( D(f ) { D(f )f}} 

^ XX X ^ X X X 

+ 6 { D ( f ) { f D(f)}} + {f D 4 (f) } + 6 {f {d ? ( f )f}} 

XX X X X 

+ 3 {f <f f 2 > } + 4 {f { D(f ) D(f)}} + {f {f D 2 (f )}} 

X XX ^ X X XX 

+ {f x { f x {y}}} J 
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x ,X - D 7 ( f ) + 21 { ( 1 )f} + 105 <l*(f )f 2 > + 35 { l/ <f ) 13(f)} 

X XX X 

+ 105 <I3(f + 210 <132(1 )f D(f)> + 35 { 13? (I ) I) 2 (f)} 

xxx ^ xx x 

+ 35 { I) 3 < f ) (f f}}* 105<f f 2 D(f )> + W <13 (f ) ( 13(f) ] 2 > 

X X xxx XX 

+ 105 <n(f )f n 2 (f)> + 105 <D(f ) f { f f» + 2 1 { D 2 ( f ) D? ( f ) } 

XX XX X X 

^ (>3 { n 2 (f ) { I3(f ) f>> + 21 { D 2 (f ) {f D(f)}} + 35 <f D(f)li(f)> 
XX XX xx 

+ 35<f xx [3(f) { f fl> + 21<f xx f D 3 (f)> + 63<f xx f { I3(f x )f}> 

+ 21<f f { f I3(f)»+ 7 {D(f ) D*(f)} + 42 {D(f ) {^(f )f}} 

+ 21 { I)(f )<f f 2 >) + 28 { D(f ){ D(f ) D(f)}} + 7 { D(f ) {f E3 2 (f )>} 

X ^ XX ^ X X XX 

+ 7 (D(f ) (f (f f»}+{f f) 5 (f)}+10{f {tf(f )f}} 

xxx" X X X 

+ 15 {f <D(f ) f2> } + 1 0 { f { I3 2 ( f ) 13(f)}} + 10 {f <f f D(f)> } 

X ^ XX XX X XX 

5 { f { 13(1 ) 1 3 2 ( i ) } } * r»{f { D(f ) {f f}}} + {l' {f 13T>(f)}} 

xx XXX X X 






(13) 
(con. ) 


+ 3 f f (f { 13 ( f )f}}} ff (f {f 13(f)}} } 

XX x xxx 

x X = 13 8 (f ) + 28 {D G (f )f} + 210 <D 4 (f )f 2 > + 56 {D 5 (f ) D(f)} 

X XX X 

f 420 <D 2 (f )f a > + 560 <D 3 (f )f D(f )> + 70 { D 4 (f ) (f )} 

-< 7 0 { D 4 ( f ) {f f}} + 105 <f f*> + 840 <D(f )f 2 D(f)> 

X X XXXX XXX 

+ 280 <D 2 (f' xx ) [D(f)] 2 > + 420 <D 2 (f xx >f D 2 (f)> + 420 <D 2 (f xx )f { f x f }> 

4- 56 {d 3 (f ) D 3 (f)} + 168 { D 3 (f ){D (f )f}} + 56 { D? (f ) {f D(f)}} 

X X X XX 

J 
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+ 280 <f f { D(f) ] 2 > + 210 <f f 2 rf(f)> + 210 <f f 2 {ff» 

^ XXX XXX XXX X 

+ 280 <D(f ) D(f) D 2 ( f )> + 280<D (f ) D(f) {f f» 

XX XX X 

+ 168 <D(f )f D 3 (f)> + 504 <D(f )f { D(f ) f » 

^ XX XX X 

+ 168 <D(f )f{f D(f)}> + 28 { ID 2 (f ) D*(f)} 

XX X X 

+ 168 { D 2 (f ) {n 2 (f )f}} + 84 {D 2 (f ) <f f 2 >) 

X X X XX 

+ 112 { i9(f ){n(f ) n(f)}} + 28 {d 2 ( f ){f D 2 (f)}} 

XX XX 

+ 28 { D 2 (f ) {f {f f})} + 56 <f DU) O 3 (f)> 

X X X XX 


168 <f D(f) { D(f )f}> + 35 <f [D 2 (f)] 2 > 

w v N vv 


+ 70 <f tftt) {f f}> + 56 <f 0(f) {f D(f)}> 

XX X N XX X 

4 35<f { f f} 2 > + 28 <f f D 4 (f )> + 168 <f f{D 2 (f)f}> 

XX X ’ XX XX X 

+ 84 <f f <f f 2 » + 1 12<f f {D(f ) D(f )}> 

^ XX XX XX X 

+ 28 <f f{f D 2 (f )}> + 28 <f f {f {f f}}> + 8 {D(f ) D 5 (f )} 

XX X XX XX X 

+ 80 {D(f ) {tfd )f}> + 120 {D(f ) <D(f )f 2 >} 

XX X XX 

+ 80 { D (f ) { D 2 (f ) D(f)}} + 80 {D (f )< f f D(f) >} 

x x X XX 

+ 40{D(f ) {D(f ) D?(f)>> + 40 { D(f ) {D(f ) {f f}}} 

XX XXX 

+ 8 { D(f ){f D 3 (f )}} + 24 {D(f ) {f {D(f )f>}> 

XX XXX 

+ 8 { D ( f ) { f {f D (f)}}} + {f D 6 (f)} + 15 {f {P*(f )f}} 

XXX X XX 



s 



+ 45 {f <D 2 (f )f 2 > } + 20 { f { D 3 (f ) D(f)}} 

x XX XX 

+ 15 {f <f f 3 > } + 60{f <D(f )f D(f)> } 
x xxx x xx 

+ 15 {f { D2(f ) tfif)}} + 15 {f {D 2 (f ) {ff}}} 

XX xxx 

+ 10 { f <f [ D( f ) ] 2 > } + 15 { f <f f D 2 (f)> > 

x XX x XX 

+ 15 {f <f f(f f}> } + 6{f { D(f ) rf(f)}} 

X XX X X X 

+ 18 { f { D(f ) { D(f )f}}} + 6{f {D(f ) {f D(f)}}} 
xxx xxx 

+ {f {f D 4 (f )}} + 6 {f {f {D z (f)f}}}+ 3 ( f x < f xx f2 > ^ 
x x xx x x x xx 

+ 4 { f {f { D(f ) D(f)}}} + {f {f {f D 2 (f)}}} 
xxx xxx 

+ {f {f {f {f f}}>> • 1 

X X X X 


(13) 

( con.) 


Corresponding expressions hold for the total derivatives of y. 

Next, we have to expand (3) in a Taylor series. Let us list the result 
for k — 1: 


fi 


- f + D(f) Cjh + ^rflf) c^h ! + | {f x f} (f) a ! ,h s 

+ | {D(f x )f} c\ h ! D«<f) c< h< + |{D*(f x )f> a\ h< 

+ 1 «*/> “> h ‘ + 715 rf(,) a ' hS + ~k { ^ H x )f} 1,5 
+ ] <D(f xx ) f 2 > a\ h 5 + ^ D‘(() «5 h' + ^ {oMOf} a« h« 

+ -h a ‘ h ‘ + 48 <W> “■ h * + ist rfm h ' 

+ ilo ( Di(t x )f> h ' + 75< tf (, xx )t! > “I 1,1 + li < D<f xxx )f)> “1 h 
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+ 01 ,f) h ’ + 7S5 f D ' (f x )f ’ “? h * < D ‘ (f ra )|S > “5 1 * 

+ -k <'* < W e > “i 1 W4 < f xxxx i< > «! h * 



Because of (4), we replaced in (14) the coefficient y by — a 2 . Using 

(4), we similarly eliminated y (u = 2, 3, . . . ) in the following expan- 

sions for ^ 2^3 • - - . We now list the expansion for f 2 . It consists of 
similar terms as (14) which are obtained from (14) replacing by a 2 
and of additional terms. We list only these additional terms: 

f 2 = •• • + {f D(f)} (y 21 a t ) h 3 + {D(f ) D(f) } a 2 (y 21 a,) h 4 
x x 

+ i {f x (y 2l "*> h * + i {f x {f x f}} (y21 “i* h4 

+ - { ° 2 (f x ) D(f)} a\ (y 21 (Vj) h 5 + D 2 (f)} a 2 (y 21 a j) h 5 

+ 2 < f x f » (r 2 i <*i> h 5 + \ D(f)> a\ (y 21 «i) h 5 

+ ~{f, D 3 (f)} (721 a 3 ) h 5 + “ {f ( D(f ) f}} (y 21 a 3 ) h 5 

u X Z X X 

S~(ir>) 

+ ~ (d 3 !^) D(f )} a 2 (y 21 a,) h 6 + ^ {^(M [^(f)} a\ (y 21 a\) h 6 f 
+ |{D 2 (f x ) {f x f» a\ <y 21 o 2 ) h 6 + |<D(f xx )fD(f)> a\ {y u a x ) h 6 
+ -^{D(f ) D 3 (f)} a 2 <y 2 i «l) h s +|{D(f ) {D(f )f}> a 2 (y 21 a 3 ,) h 6 

DA Z X X 

+ 4 ^ f xx f 012 (y 2i «l) h6 + \ < f xx f { f x f }> a 2 ( T21 «i) h e 

+ 2 <f xx [D(f)]2 > <?2l *i) 2 h 6 {f x D*(f)> (721 «i> h 6 
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4 ^ {f <f f®> } <y 21 n\) h 6 + \ {f { D 2 (f )f}] (>'2i A) h 6 
8 x xx 4 x x 

IXO) 0t\ (y 21 rv^hU-^lD 3 (O D 2 (i)> n\ (y 2( a \) h? 


7^ {^( f x ){ f x f » <4 <y 2 , a' 2 1 )h 7 + |<D 2 (f xx )f D(0> <4 (y 21 a,) h 7 


-7- {l^(f ) n^(f)} a 2 (y 21 a 3 )h 7 + j { D 2 (f ) f D(f ) f} } a 2 (y 21 a 3 ) h 7 

1 Z X 4 X X 


+ 4 < ^ D(f XX )f ° 2(f) > a 2 (y 21 «l)h T+ | <«( f xx )f {yl >02 (yfel O 2 )h 7 


T <D(f ) [D(f)] 2 >a 2 (y 21 o 1 ) 2 h 7 + 7<f f 2 D(f)> « 4 2 (y 21 aj)h 7 

Z XX ” XXX 


8 ^ xxx 


2^ D(f x ) n i ^2! o 4 )h 7 + ^{D(f x ) <f vv f 2 >) »2 (yzi^l) h 


8 x xx 


4{o(fJ { D 2 ( ) f } } « 2 (y 21 o 4 1 )h 7 + ^<f xx f u 3 (f)> cv\ (v 21 a 3 )h 7 


X X 


4 7<f f { r>a )f)> rvl (y 21 ai)h 7 + 7<f D(f) D 2 (f)> (y 21 afi) ( Y21 <* 2 ) h’ 

4 XX X 6 XA 


+ |<f xx D(f) {y}> (y 21 a,) (y 2 io 1 )h 7 +T^{f x D 5 (f)} (y 21 c 4 )h 7 


+ |{f <o(f )f 2 » (y 21 o 5 !)h 7 + -^(f (rf(f v )f)} (y 21 «i) h 7 

n X XX 1 Z X X 


+ ~~ { D 5 (fy D(f ) } af (y 21 ai)h 8 + {dMO D 2 ^)} a 4 (y 21 at 2 ) h 8 


+ ^ { D‘ ( f x ) ( f x f}} o 4 2 <V2i * 2 i)h 8 + <rf(f xx )f D(0> (y 21 0!)h 8 


+ — (d* (f ) D 3 ( 0 ) a 2 (>« a’) h 8 + ~ { D? (f ) { D(f )f }} a\ (y 21 a 3 j) h 8 
36 x x x 
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A 


+ 8 < t > 2 < f xx )f 02 <721 + 8 ^ I:)Z<f xx )f ^ f x f} > a * (yu Q, ‘ )h8 

+ “ <D 2 ( f )[D(f)] 2 >rv| (y 2 l «l) 2h8 + i < D <f ) f2 D(f )> «2 <721 “l> h8 

4 XX O XXX 

+ 7r { D 2 (f ) D*(f)> Oi\ <721 «D h' + ^{rf(f ) {tf(f )f}} ot\ <721 Q 'i)h 8 
4 O X O X X 

+ -77: { D 2 ( f ) <f ,f 2 >) c y\ <721 «1)h 8 + -k <D(f )f D»<f)> a\ ( 721 ^)^ 

It) X XX lei XX 

+ 4 < D(f xx )f{D(f x K}> "* <T2l « 3 i)h 8 
+ ^ < D < f xx > D < f ) ^<f)> “2 <721 Q l) <721 a\) h 8 

+ | < C D < f xx ^ D < f ) {y» a 2 <721 «i) <721 <**) h ° 

+ ~lfi < f xxx f2 ° 2(f) > "2 <V 2 i « 2 ) h* 

+ ^ < f xxx f2 {f x f}> ° 2 (T2 ‘ ft2) ^ + 4 <W 1 U(f)]2 > «2 <V2i «l) 2 h 8 
+ "2 <7 2 i «5)h 8 + ^{D< f x ) {D 3 (f x )f}} a 2 (7 2 ia 5 i)h 8 

+ 8 ^ D(f x ) < D<f xx )f2 » " 2 (T21 CV ‘ )h8 + 48 < f XX f D<(f) ^ > “2 <721 a 4 i>h 8 

+ 8 <f xx f f° 2(f x )f y a 2 <721 «l)h 8 + <f xx f < f xx f2 ^> a \ <721 "l) 118 

+ ^< f xx D < f ) D ? <f)> <721 a,) <721 a 3 i)h 8 

+ |< f xx D «> (D<f x )f}y (721 «i> <721 «i) h” 


y (i5) 

C (con.) 


J 
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+ -i<f [D ! (f)] 2 > <Yn«l> V+ 7<< rf(f){y» (YH“ ! .) !h * 

o XX 


+ ^<f {f f} 2 > (y 2 i« 2 i) 2 h 8 + ^{fD 6 (f)} (Y 21 « 6 i)h 8 

8 xx x ‘■IV x 

+ -Mf { D* (f )f}} (72iffl)h 8 +77{f <&(f )f> 

48 x x 10 x aa 

+ Ii {f K <'**/>> ( ^" 5)h> 


(15) 

( con. ) 


Next, we list the expansion of f 3 . This expansion consists of terms simi- 
lar to those of (14); these terms are obtained replacing a 1 by o 3 in (14). 
Furthermore, f 3 contains terms similar to those of (15); these terms are 

obtained replacing in ( 15) ot\ by a% and (y 21 b y O 31 a t + 732 a l )• 
Finally, f 3 contains additional terms which do not appear in (14) or (15). 
We now list these additional terms of f 3 : 


. + {f (r DU))) 732 <7*1 <Vl)h 5 + ( D<1 ){f ixtm «S • 7» (721 a h * 

X X 

< If { n(f ) D(f)}}7ii«* (72l«l) h6+ i{ f x {f x Il2U)}}y32 (y2 * ft > )h6 

x X 

+ i {l x l,(n]) " 3 ‘ T32 (y2 ‘“‘ )h7 

+ (nu ) { 1 )( 1 ) D(i)}} «j • Y32 (721 

X x 

>( 16 ) 

+ -{[)(f ) {f I J 2 1 1 > } > n 3 • 732 (72J “ 2 l) h? 

2 x x 

\ - { l)(f > { f ( I' f)}} «3 • 732 (721 h? 

2 x xx 

-1 - <f f {f I > ( f ) » 03 ' 732 <721 ^l)* 17 
2 xx x 

+ i{f <D 2 (l ) D(i)}}73 2 «I(72.«l)h 7 +T{l x <D(f K ) D 2 (l'))}(732«2) h? 

2 x x 


J 
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{ D(f ) {f f}}} yi2 <X 2 (?21 «l)h 7 + | U <f f D(f)> ) (732 «i)^21 

Z X X X ^ A AA 


+ ^{f {f D 3 {£) }} r 3 2 (Y 21 «i)h 7 + |{f {f {D(f x )f}}}r32 (y 2 i«i)h 7 

6 X X Z X X A 

+ -{£ { D? (f ) D(f )}} (y 32 «2> ^a^h 8 +|{f ( {D 2 < f x ) D 2 < f )}} <Y 32 a|) (y2i«i)h 8 
6 x x 

+ 7 {f { D 2 (f ){f f}}} <y 32 a\) (y 2 i“i)h 8 + |{f <D(f )f D(f)> } <y 32 a?) <y 2 i a^h 8 

4 X X x Z X XX 

+ 7 {f { D(f ) tfU)}} (y 32 a 2 ) (y 21 a 3 )h 8 + | {f (D(f ) { D(f )f}}> (y 32 a 2 ) <Y 2 i a 3 i)h 8 

6 X X ji X X X 

+ | {f x < f xx f D 2 (f)> } (y 32 a§) (y 21 a^h 8 + ~ {f x <f xx f (f x f }> } (y 32 ajj) (y 2 i a 2 ,)h 8 
+ | {f x < f xx I d( 0] 2 > }Y32 <y 2 i« 1 ) 2 h 8 + ^{f x {f x D‘(f)}}y 3J (y 31 a*)h 8 

+ |{f x (f x {D 2 <f x )f}»yj 2 <Y 2 i a 4 t )h 8 + 1 {f x {f x <f Qt ^> }} y 33 (y 2 i «1>h 8 
+ <f D(f) {f D(f)}> y 32 (y 3 , + y 3 2 «j) (y 21 ai)h* 

^ XX X 

+ — <^f xx f { D(f x ) D(f)}^> a 3 (y 32 a 2 ) (Y 21 a l^ 

+ |< f xx f { f x rf(0}> ®! • Y32 (Y2 1 a 2 1 )h 8 +|<f xx f{f x {f x f }» a 2 Y 3 2 ( Y21 a 3 ) h 6 

+ |{D<f x ) {rf(f x ) D(f )}} a 3 (y 32 a 2 ) (y 3 i a 3 )h s 
+ |{D(n{D(f x ) tf(f)}} a 3 (y 32 a 2 ) <Y 2 i a 2 j)h 8 

+ ^ { D(f x ) { D(f x > {f x f}}} a 3 (y 3 2 a 2 ) (Y 21 a 2 )h 8 

+ |{ D(f x > <f xx f D<f)> } a 3 (y 32 a\) (y 21 a,)h 8 

+ “ { D(f ) {f D^f)}} a 3 • y 32 (y 2 i a 3 j)h 8 +|{D(f ){f {D(f )f}}} a 3 • y 32 (y 21 a^h 8 
6 x x z x x a 

+ 1 {d 2 ^) {D(f x > D(f)}} a 3 (y 32 a 2 ) (y 2 i a,)h 8 + ^ (D 2 ^) {f x O 2 (f >1} a 3 • y 32 (y 2 i a 2 )h 8 


(16) 

( con. ) 
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+ — { D 2 (f ) {f {f f}}} «3 • Y32 <T 21 «l ) h8 

4 . x x x 

+ -<D(f )f{f D(f)}> • y 32 (y 21 ai)h 8 

2 ^ XX X 

+ -{D?(f){f D(f)}} al ■ y 32 (y 2 i ai)h 8 
6 xx 

Finally, we list the additional terms for f 4 : 

f. = ...+{f {f {f D(f )}}} y 43 • 732 O 21 “i ) 1 * 7 \ 

4 X X X 1 

+ {f { D(f ){f D(f)}}} (y 43 «s> Y32 <y 2 i “l)* 1 ® I 

1 X x X I 

+ {f {f { D(f ) D(f )}}} y 43 (y 3 2 “ 2 ) <Y 2 i«i > h8 I 

XXX \ (17) 

+ -{f {f {f D^lf)}}} T 43 * "F32 ^21 «'l> h8 I 

2 x x x I 

+ -{f {f {f {£ f}}}> 743 • 732 (721 “l > h8 I 

2 x x x x I 

8 / 

+ { D ( f v ) { f {f D(f )}}} a 4 y 43 732 (721 «l) h / 

1 x X X / 

When proceeding to f 5 , ... , no more additional new terms for h 8 or 
lower powers of h are obtained. 

Formulas corresponding to (14), (15), (16), and (17) hold for gj, g 2 , 
g3) g4 . in all these formulas, all functions on the right-hand side are 

to be taken for t = t 0 . 

7 Introducing the expressions for f t , f 2 , - • • i nto (5) and comparing the 
resulting expansion (5) term by term with the expansion (8), using the 
values (13) for the total derivatives in (8), we obtain the equations o 
condition for the Runge-Kutta-Nystrom coefficients as listed in Table 1. 



1. All tables are at the end of this report. 
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The left-hand part of Table 1 represents the equations of condition for 

x, y as obtained from the expansions of (3) and the first (or first two) 

equations (5); the right-hand part of Table 1 states the equations of 

condition for x, y, as obtained from the expansion of (3) and the last 

equation (5). For the right-hand part of Table 1, the weight factors 

c of the table are to be replaced by the weight factors c of the first 
k k 

derivatives x, y. In the marginal columns of Table 1, one finds listed 

the h-power of the corresponding term in the Taylor expansions for x, y 

(left-hand side) or x, y (right-hand side). These terms are listed in 

the same order as in equations (13). Some of the equations of condition 

of Table 1 are duplicates or multiples of other ones and are marked by an 

asterisk (*). When we consider the leading term of the local truncation 

error, these duplicates must also be taken into account. 

Later in this paper, we will refer to the equations of Table 1 by a Roman 
numeral (I through X), indicating the h-power (left or right marginal 
column in Table 1) to which the equation belongs and by an Arabic 
numeral indicating its position in the block of the equations of the h-power 
under consideration. For example, the last equation in Table 1 should be 
referred to as (X, 72) or, if considered as an equation with the weight 
factors c , as (IX, 72)* . 

K 

Equations (X, 1) through (X, 72) of Table 1 are required only as equa- 
tions (IX, 1)" through (IX, 72)* for the local truncation error in x, y. 


8. We reduce the equations of condition of Table 1 by omitting the duplicate 
or multiple equations ( * ) and by introducing, in addition to (6) and (7), 
the following assumptions: 



16 



Y21 “i = jr 2 "2 


731 “l + 732 a l - 12 "3 

741 “l + 742 "I + 743 "3 = ^2 ^ 

7in a l + 7112 «2 + 7113 ^ I + + 7ino “ io =^2 ““ 

741 «1 + 742 “1 + 743 "3 = 20 
7 5 1 “l + 752 “I + 753 “I + 754 a l = ^ “I 

\ 

/ 

Till « 3 i + 7ii2 «2 + 7113 “i + • • • + 7ino “10= ^ “11 

7ioi 

c 4 741 + c 5 y 51 + . . . + c 9 y 91 + c 10 < = 0 

7ui 

c 4 74i + ^5 75 1 + • • • + C 9 y 91 + c 10 y 101 = O 

L 7 101 I 

c 4 «4 y 4 j + C 5 a 5 7 5 1 + . . . 4- C 9 a 9 y 91 + c 10 i ( = ° 

7m | 

c 4 a A y 4 j + c 5 a 5 y 5 i + . . • + c 9 a 9 y 91 + c 10 y 101 = 0 

c 4 a\ 741 + c 5 “1 751 + * - * + c s “I 7gi + C 10 7m = 0 
c 4 a\ y 4l 4- c 5 a 2 5 y 51 + . . . + c 9 a\ y 91 + c 10 7ioi = 0 j 



( 20 ) 


( 21 ) 


( 22 ) 


( 23 ) 


( 24 ) 
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(25) 


c 4 742 + c 5 752 + • . . + c 3 y 92 + C 10 y u2 - 0 I 

^4 T 42 + ^5 752 "*'••• + ^9 "K 92 + C 10 Ti 02 = 6 j 


c 4 743 + c 5 Y 53 + • • • + <-'9 y fl3 + Cjq y 113 - 0 
^4 743 + ^5 T53 + • • • + ^9 T93 + ^10 7 103 = 0 


(26) 


e 5 V 5 4 741 + c 6 <Tfi4 741 + res 751 ) + 


t: 5 754 741 + ^6 <764 741 + 7fi5 7s 1 ) + 


+ c 9 (y 94 y 41 + . . 
+ c 10 ^7114 741 + • 
+ c 9 (794 741 + • • 
+ c 10 (y 104 741 + ■ 


• '■ 798 78 1) 

• • + 7 1110 7 ioi) = 0 
- + 7»8 781 ) 

• • + 7 109 791) - ^ J ■ 


> (27) 


In the first equations (22) and (23) , the upper line (y 101 ) holds for the 
eighth-order formula and the lower line ( y UJ ) for the ninth-order formula. 
Introducing these assumptions, we convert the necessary and sufficient equa- 
tions of condition of Table 1 into a system of sufficient equations of condition 
that can be solved in a relatively easy manner. In fact, these assumptions 
transform the equations of Table 1 into several separate systems of linear 

equations for the coefficients y 

kX 


Let us now consider the reductions of the equations of Table 1 effected by 
the above assumptions. Assumptions ( 19) — together with the first 
equation (18) — lead to the following identities in Table 1: 


(V, 3) s (V, 1), (VI, 4) = (VI, 1), (VII, 4 ) = (VII, 1) 
(VIII, 4) = (VIII, 1), (VIII, 9) s (VIII, 1), (IX, 4) = (IX, 
(IX, 10)3 (K, 1), (X, 4) = (X, 1), (X, 11) E (X, 1) 

(X, 35)3 (X, 16), (X, 61) 3 (X, 56) . 





(28) 


J 


Because of the identities (28) the equations (V, 3), (VI, 4), (VII, 4), 

(VIII, 4), (VIII, 9), (IX, 4), (IX, 10), (X, 4), (X, 11), (X,35)’, (X.fil), 
and the corresponding equations in c can be omitted from Table 1. 

K 
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Assumptions (20) — together with (19) and the first equation (18) — yield 
the following identities in Table 1: 


(VI, 5) = (VI, 1), (VII, (i ) = (VII, 1), (VIII, 7) = (VIII, 1),^ 

(VIII, 19) = ( VIII, 18), (IX, 7) = (IX, 1), (IX, 1(5) = (EX, 1) 

(IX, 25) h (IX, 24), (IX, 32) - (IX, ,20), (X, 7) = (X, 1) , V (29) 


(X, 20) = (X, 1), (X, 29)= (X, 28), (X, 22) s (X, 1) 

(X, 48)= (X, 47), (X, 59)= (X, 58), (X, 71)= (X, 70) , y 


eliminating equations (VI, 5), (VII, 6), (VIII, 7), (VIII, 19), (IX, 7), 
(IX, 16), (IX, 25), (LX, 32), (X, 7), (X, 20), (X, 29), (X, 33), 

(X, 48), (X, 59), and (X, 71) from Table 1. 


Assumptions (21) , together with (19) and the first equation (18), lead to 
the following identities: 


(VII, 8) = (VII, 1), (VIII, 12) = (VIII, 1), (IX, 13) = (IX, 1),[ 
(X, 14) 5 (X, 1), (X, 31) = (X, 1), | 


(30) 


thereby eliminating equations (VII, 8), (VIII, 12), (IX, 13), (X, 14) 
and (X, 31) from Table 1. 


Because of ( 7) and the first equations (18) and (22) : 
Till = Yioi = 0 

Similarly because of (7) and (18) : 


Yil2 - °> Yii3 - 0 • 


(32) 
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Therefore, we can omit the last term on the left-hand side of (22) , (23) , 
and (24) and of the first equations (25) and (2(5). We can also omit the 
last term on the left-hand side of the first equation (27) since (7) and the 
first equation (22) hold. 


From (22) and the previous assumptions, the following identities result: 


(VII, 10) = (VII, 8), (VIII, 18) = (VIII, 10), 
(IX, 30) s (IX, 27), (X, 56) s (X, 53) , 


(33) 


eliminating equations (VII, 10), (VIII, 18), (IX, 30) , and (X, 56), from 
Table 1. 


Assumption (23) and previous assumptions lead to: 

(VIII, 14) = (VIII, 12), (IX, 24) = (IX, 21), (X, 46) s (X, 43) , (34) 

eliminating equations (VIII, 14), (IX, 24), and (X, 46) from Table 1. 

By (24) and previous assumptions equations (IX, 15) and (X, 28) are 
eliminated from Table 1 since they then become identical with (IX, 1) 
and (X, 25), respectively. 

Finally, assumptions (25) and (26), together with previous assumptions, 
lead to the identity: 

(IX, 34) = (IX, 27), (X, 64) = (X, 53), (35) 

and assumption (27), together with previous assumptions, to: 

(IX, 36) = (IX, 27), (X, 66) = (X, 53), (X, 70) = (X, 67), (36) 

thus eliminating equations (IX, 34) (IX, 36), (X, 64), (X, 66), and 
(X, 70) from Table 1. 

Omitting in Table 1 all equations which are duplicates or multiples of other 
ones or which can be eliminated by using the above assumptions, the table 
reduces to the following equations: (II, 1), (III, 1), (IV, 1), (V, 1), (VI, 1), 
(VII, 1), (VIII, 1), (IX, 1), (vm, 15), (IX, 21), (IX, 27), which are 
listed in this order in Table 2. To write the last three equations in a more 
concise form, we introduced in Table 2 the abbreviation: 
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(37) 


\i a i +y M2 *2 


+ y , cat = P 

/i,J*-l M-l 


The remaining equations in Table 1 which belong to ninth-order terms for 
x will be considered later when dealing with the local truncation errors in 
x and x. 


In the ninth equation of Table 2, the upper line (P 104 ) holds for the eighth- 

order formula and the lower line (^-) for the ninth-order formula. From 

o u 

the two lines of this equation, it follows that 

_ J_ 

PlM ~ 30 


must hold. 

The equations of Table 2 have to be solved together with the assumptions 
( 19) through (27) . 

9. The first eight equations of Table 2 do not contain the coefficients y . 
They represent linear equations for the weight factors c^ or c^. By 
solving these equations, we can express the weight factors by the a^'s. 

Next, we have to solve for the coefficients y . Obviously, equations 

K A 

(22), (23), (24), and (27) are satisfied by 

T 41 = T51 = Tel = 771 = V 8 i = Y91 = T101 = Till = ^ 

The first equation (19) together with the first equation (20) leads to a 
restrictive condition for oq: 


Similarly, the third equation ( 19) , together with the third equation (20) 
and the first equation (21) , yields because of y 41 = 0, the following re- 
strictive condition for a 2 : 



5a a ~ 3oq 
2a 3 - a A 


(40) 
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The remaining a 's may be chosen arbitrarily. We tried to select these 

K 

-values in such a way that the leading term of the truncation error be- 
comes small for our Runge-Kutta-Nystrom formula. 

The remaining coefficients y can easily be determined as follows: The 

K A 

first equation (19) yields y 21 ; the second equation (19) and the second 
equation (20) determine y 31 andy 32 . Since y 41 - 0, from the third equation 
(19) and the third equation (20) the coefficients y 42 and y 43 can be obtained. 
Because of y 5 j - 0 we can obtain the coefficients y 32 , Y53> and y 34 from the 
fourth equation (19), the fourth equation (20), and the second equation 
( 21 ). 

f rom the ninth and the tenth equations of Table 2, we can determine Pg 4 , 
P74 > P 84 » anc * P94 * since these equations are four linear equations for 
these quantities ^ ^44 ’ P54 being known already). 

The fifth equation (ID), the fifth equation (20), the third equation (21), 
and P B4 , together with y 61 = 0, yield the coefficients y 62 , y 63 , y 64 , and y B5 . 

Since we have more coefficients y than equations of condition, we may 

K A 

set some of these coefficients equal to zero: 

772 = 782 = T83 “Tim = (41) 

The coefficients y 73 , y 74 , y 75 , y 76 can then be obtained from the sixth 
equation (19), the sixth equation (20), the lourth equation (21), and P 74 . 
Similarly, y 84 , y 85 , y 86 , and y 87 are determined by the seventh equation 

(19) , the seventh equation (20), the fiftli equation (21), and P 84 . 

Next, we use the first equation (25) and the first equation (26) to find 
and y 93 . The first line of the last equation in Table 2 yields P_, 5 . This 
value P 95 , together with P 94 , the eighth equation (19), the eighth equation 

(20) , and the sixth equation (21), determines the coefficients y^, 795, y 98 , 
y 97 , andy 98 . Similarly, we use the second equations (25) and (26) to 
compute y 102 andy 103 . From the second line of the last equation in Table 

2 we find P 105 . Since y 104 = 0, we can determine the remaining coefficients 

7105* Tioe, 7 107 , 7i08> and y 109 from P 105 , P 10( (-— ), the ninth equation (19), 
the ninth equation (20), and the seventh equation (21). 
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This concludes the computation of the coefficients y ^ since y 

(A= 0, 1, 2, 10) is given by (7), and y (k = 1, 2, 11) can 

K U 

be determined from (4). 

10. Table 3 lists the coefficients of an eighth-order Runge-Kutta-Nystrom 
formula RKN 8(9) . The expression TE in Table 3 represents the 

X 

leading term of the local truncation error in x for our eighth-order 
formula and is obtained as difference between the first and the second 

formula ( 5) . 


In Table 4, the error coefficients for the leading truncation error term 
and x and x are listed. If the error coefficients differ by a constant 
numerical factor only, the coefficient with the largest factor is listed. 
The product of the error coefficient and the corresponding expression 
in the partial derivatives, as listed in (13), times h 9 is the actual con- 
tribution to the leading truncation error term. 

Because of P lft j = — , the tenth equation of Table 2 holds also for the 
1U4 30 

eighth-order formula RKN 8, so that this equation does not contribute to 
the leading truncation error term in x. The only contributing equation 
of Table 2 is the last one [equation (IX, 27) of Table l] . Replacing 
equation (IX, 27) by equation (IX, 29), because of the larger factor of 
the latter one, the error coefficient that determines the contribution to 
the leading truncation error term in x reads: 


L 29 


77 <C 4 P 45 + c 5 P 


55 


+ Cc 


P95 + c 10 P 


105 


) - 


15 


362 880 


(42) 


In the case of x, there are seven essentially different error coefficients 
that contribute to the leading truncation error term: T 10 , T 26 , T 4S , T 5l , 

T 52 > ^58 > an d Tg8- 

It is essential that the error coefficients in x are of about the same order 
of magnitude as the error coefficients in x. If the error coefficients in 
x are large compared with the error coefficients in x, large truncation 
errors in x might be generated since in our method the integration step - 
size is determined by the local truncation error in x (by TE^) • 
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SECTION 1 1. SEVENTH-ORDER FORMULA RKN 7(8) 


11. The derivation of a seventh-order Runge-Kutta-Nystrom formula RKN 7(8) 
closely follows the same pattern as in the case of our eighth-order formula 
RKN 8(9) of Section I. Naturally, a seventh-order formula requires fewer 
evaluations of the differential equations per step than an eighth-order for- 
mula. We shall present in the following a seventh -order formula based on 
nine evaluations per step, a tenth evaluation being taken over as first eval- 
uation for the next step. 


In the case of a seventh-order formula, equations (3), (5), and (6) must 
then be replaced by: 

f 0 = f(t 0 , x 0 , y 0 ) 


K- 1 


t = f xo + xoo^h + h 2 . 


A=0 


K-l 

y„ + y, ah + ll ! • V y Kk g x 

A— U 


8 


(k = 1, 2, , 9) 


X = x 0 + x 0 


+ Xnh + h 2 ■ y c f + 0 ( h 8 ) 

n K K 


9 

= x 0 + x 0 h + h 2 • V £ f + 0 ( h 9 ) 
w u j K K 

K = 0 


(43) 


(44) 


8 

x — Xq + h • y c f + 0(h 8 ) 

n K K 
K = 0 


c = c for k = 0, 1, 2, . . . , 7 

K K 


c 8 = 0 


(45) 


A 

c 9 = c 8 
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In Table 1 the last 72 equations of condition now have to be omitted, 
since they correspond to tenth-order terms in the Taylor expansion for 
x or to ninth-order terms in the expansion for x. 

Instead of equations (18) through (27), we now make the following 
assumptions! 


A A . . 

Cl = Cj = c, = 0, c 2 = c 2 = c 2 = 0, a g = or, - 1. 


(46) 


1 3 

721 a l ~ g rt 2 


731 °T + 732 a 2 ~ g a 3 


1 3 

791 Q 1 + 792 ^2 + • • • + 798 a 8 " 7 °9 


A 


> 


(47) 


731 "1 + 732 a 2 


2 _ — Qfi 
" 12 3 


A 


741 Q'l + 742 Q'l + 743 - 12 Q* 




791 Q Z 1 + 792 a 2 + • • • + 798 a S - 12 "9 J 


7*i Q'l + 732 a l = ^ Q3 5 


741 Ql + 742 a 2 + 743 a 2 ~ 20 






791 Q 3 + 792 “2 + . • • + 798 Qs 


20 


at 




(48) 


(49) 
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c 3 731 + c 4 741 + • • • + c 7 771 + c 8 

C3 731 + c 4 Y41 + • • • + c 7 771 + c 8 y 81 = 0 J 

C 3 ^3 731 + c 4 a i 741 + • • • + c 7 «7 771 + c 8 7 9 1 =■ 0 

c 3 «3 731 + C4 «4 741 ■*■••• + c 7 «7 771 + c 8 781 = 0 

c 3 732 + c 4 742 + • • • + c 7 772 + c 8 792 = 0 

C3 732 + C 4 742 + • • • + C7 772 + C 8 782 “ 0 j 

Because of y 91 = Cj = 0, we obtain from the first equation (50): 




(50) 


(51) 


( 52) 


7 8 i=0 . (53) 

The assumptions (46) through (52) reduce the equations of conditions 
of Table 1 to those of Table 5. 


12. We now have to solve the equations of Table 5 together with the assump- 
tions (47) through (52): The first seven equations of Table 5 allow us 
to express the weight factors by the a's. From the first equation (47) 
we obtain y 21 . The second equation (47), the first equation (48), and 
the first equation (49) lead to a restrictive condition for the ads: 


O' 



So 1 ? - 302 

2a 2 - Oj 


(54) 


and to 7 3 j, 732> The third equation (47), the second equation (48), and 
the second equation (49) yield y 41 , y 42 , y 43 . 

On the other hand, equations (50) and (51) allow us to express 741 , y 51 , 
Y 8 l, and 7 ? i by 731, the cv's and the weight factors. Therefore, we have 
two different expressions for 741. Equating them leads to another 
restrictive condition for the cv's: 
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(2« 3 - a 4 ) 


(55) 



Since y 51 is already known, the coefficients y 52 , y 53 , Y 54 are obtained 
from the fourth equation (47), the third equation (48), and the third 
equation (49). 

Putting y 62 = 0,we obtain, since y B1 is already known, y 63 , y M , Yes from 
the fifth equation (47), the fourth equation (48), and the fourth equation 

(49). 

Since Y 92 = c 2 = 0, the first equation (52) yields y 72 . Since y 71 is known, 
we can obtain y 73 , y 74 , Y 7 s> 776 from the sixth equation (47) , the fifth 
equation (48) , the fifth equation (49), and the last equation of Table 5 
(written as equation with weight factors c^). 

Similarly, the second equation (52) yields Ys2* Setting y 83 =0, we obtain 
y M , y 85> y 8g , y 87 from the seventh equation (47), the sixth equation (48), 
the' sixth equation (49) , and the last equation of Table 5 (written as 
equation with weight factors c^). 

This concludes the computation of the coefficients y^ ^ since y 

(A = 0, 1, . . . , 8) are equal to the weight factors c^ (A= 0, 1, . . . , 8) 

and y (k = 1, 2, ..., 9) can be determined from (4). 

K 0 

13. Table 6 lists the coefficients of a seventh-order Runge-Kutta-Nystrom 
formula RKN 7(8). The restrictive condition (55) makes it somewhat 
harder than in the case of the formula RKN 8 (9) to find reasonably 
simple a -values. This explains the somewhat unwieldy coefficients 
y of Table 6. 

y kX 

We notice in Table 6 that a 2 and c 4 are negative. It is possible to obtain 
positive values for all a's and all weight factors by a different choice 

of a z (e. g. , a 3 = j-} . However, the coefficients y^ then turn out to 
be even more formidable. 

In Table 7 the error coefficients of the formula of Table 6 are listed in 
the same way as we have listed in Table 4 the error coefficients of the 
formula of Table 3. 
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SECTION III. SIXTH-ORDER FORMULA RKN 6(7) 


14. For the derivation of a sixth -order It unge-Kutta-Ny strom formula 
RKN 6(7), we proceed in a quite similar way as in Sections I and II. 

In the case of a sixth -order formula RKN 6(7) we have to consider only 
the first 43 equations of Table 1. We base our sixth-order formulas on 
seven evaluations of the differential equations per step and allow for an 
eighth evaluation that will be taken over as first evaluation for the next 
step. We then have: 


f 0 = f (t 0 , x 0 , y 0 ) 


K — 1 


f K = f \ t() + a K h ’ x o + x o h + h 2 • y; y KX f x , 


y» + y«h + h ! -V y » 


X- 1 

-V 

LJ 

A-0 


(k = 1, 2, , 7) 


x - x 0 + x 0 h + h 2 • 7 c f + 0 (h 7 ) 

1 u K K 

K~ 0 


and 


x = x 0 + x 0 h + h 2 • Y c f + 0 { h 8 ) 

u -i K K 

K = 0 


X = Xn + 


h • t 6 f + 0 (h 7 ) 


K - 0 


K K 


and 


c = c for k~ 0, 1, 2, . . . , 5 

K K 


c 6 = 0 


c 7 - c 6 


(56) 


(57) 


(58) 
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In the case of a sixth-order formula IIKN 6(7), we make the following 
assumptions: 


c i = Cj = Cj = 0, a 6 = q 7 = 1. 

1 3 

721 a \ = g a i 


A 


i 

731 ff l + 732 a 2~ Q 


a l 


> 


( 59 ) 


(60) 


1 3 

771«l+ 772 0'2 + ••• + 776 a 6 = g *7 


2 1 4 

721 «1 ^ 75 a 2 


J 


731 «1 + 732 a 2 = 7^ «3 


(61) 


771 «1 + 772 ^2 + 


776 


CV 


2 _ 


12 






C 2 721 + c 3 731 + c 4 741 + c 5 751 + c 6 771 " 0 
C 2 721 + d 3 731 + 6 i 741 + c 5 y 51 + C 6 y 61 = 0 


(62) 


Obviously, these assumptions reduce the equations of Table 1 to those of 
Table 8. 

15. We now solve the equations of Table 8 together with the assumptions 

(60) , (61) , and (62). The first six equations of Table 8 yield the weight 
factors as functions of the a's. From the first equations (60) and (61) , 
we obtain the restrictive condition: 
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and the coefficient y 21 . The second equations (60) and (61) yield y 31 and 
y 32 . Setting y 41 = 0, we obtain y 42 , T 43 from the third equations (60) and 
(Gl). Since y 71 = c t --- 0, the first equation (62) yield y 51 . The coef- 
ficients y 52 , y 53 , y 54 are then obtained from the fourth equations (60), 
(Gl) and the last equation of Table 8, written as equation with weight 
factors c . 


The coefficient y 61 is determined by the second equation (62). Setting 
y fi2 - 0, the coefficients y fi3 , y c4 , y 65 can be computed from the fifth 
equations (60) and (61) and from the last equation of Table 8, written as 
equation with weight factors c . 


This concludes the computation ot the coetlicients y^ x since y^ 

(A - 0, 1, 6) are equal to the weight factors c (A - 0, 1, .... 6) 

andy (« - 1, 2, ... 7) can be determined from (4). 
k0 

1G. in Table 9 the coefficients of a sixth-order Runge-Kutta-Ny strom formula 
RKN 6(7) are presented. In Table 10 the error coefficients for the 
leading truncation error term in x and x arc listed. Again, if the enor 
coefficients differ by a constant numerical factor only, the coefficient 
with the largest factor is listed. As Table 10 shows, the error coef- 
ficients in x are of about the same order of magnitude as the error 
coefficient in x. 

SECTION IV. FIFTH-ORDER FORMULA RKN 5(6) 


17. To derive a fifth-order Runge-Kutta-Nystrom formula RKN 5(6), we pro- 
ieed in a quite similar way as in Sections I to III. We base out fifth 
order formula on six evaluations of the differential equations per step 
and allow for a seventh evaluation that will be taken over as first evalua- 
tion for the next step. We then have. 


f 0 f(t 0 ,x 0 ,yo) 


f 

K 


K-i 


f(to + ah, x 0 + x 0 <vh +h 2 - V y ( , A f A , 


A-0 

K — 1 


yo 


+ y 0 «,h+h 2 .^ y KX g x 


A=0 


(k = 1 , 2 , ... , 6 ) 



(64) 
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x = x„ + k„h + h J -£ c Mort 


x.+ k„h+h ! -£ ^ f K + 0th!) 
K^O 


k o + h -Z fc A + 0(h6) 


c — c for k - 0 , 1 , 2 , 3 , 4 

K k 


c 6 = c 5 


We make the following assumptions: 


Ci - 0 , « 5 = a 6 = 1 * 


Ci = Ci - C t = U , “5 

721 “l = 

1 3 

731 tt l + '1'32°'2 ~ ~n “3 


741®! + 742 a 2 + T43 a 3 ~ 6 ^4 


y u a i+ 752^2 + 753 “3 +754 “4 " fi a 5 

_L 3 

'y 61 0f 1 + 762 a 2 + 763 a 3+ 764 <*4 + ^65 “ 5 ' 6 


7 31 = 0 , 741 = 0 > V51 = 0 ’ 'Vet - 0 > 

which reduce the equations of condition of Table 1 to the equations of 
Table 11. 
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18. The first five equations of Table 11 yield the weight factors c or c as 

functions of the a ! s. 

K 

From the first two equations (68) we obtain y 21 and y 32 . From the third 
equation (68) and from the last equation of Table 11, this equation 
written with weight factors c , the coefficients y 42 and y 43 can be deter- 
mined. Using the fourth equation (68) and the last equation of Table 11, 
written as equation with weight factors c , the coefficients y 53 and y 54 can 

K 

be expressed by the cv ! s and by y 52 . By still having y 52 available, we 

K 

could select y 52 in such a way that some of the error coefficients of our 
fifth-order Runge-Kutta-Nystrom formula become small. By doing so, 
we obtained error coefficients, some of which were considerably smaller 
than those obtained from setting y 52 = 0. The coefficients y of the 

last equation (68) are naturally equal to the weight factors c , since it 

A 

was intended to take over the seventh evaluation as first evaluation for 
the next step. 

19. In Table 12 are listed the coefficients of a fifth-order Runge-Kutta- 
Nystrom formula BKN 5(6) , and in Table 12 the error coefficients in x 

as well as in x for the formula of Table 12. Table 12 shows that the error 
coefficients in x are of the same order of magnitude or smaller than 
the error coefficient in x. 


SECTION V. FOURTH-ORDER FORMULA RKN 4(5) 


20. A fourth-order R unge-Kutta-Nystrdm formula RKN 4(5) can be based on 
four evaluations per step if we allow for a fifth evaluation that will be 
taken over as first evaluation for the next step: 


fo f(to, x 0 , yo) 

/ K'-l 

f = f ( t 0 + a h, x 0 + xot* h + h 2 • V y f , 
K \ K K *■ — ! bC A K 


yo 


yo a h 


K — 1 

•E 

A.-0 


T kA S n 


(k 1,2, 3, 4) 


( 70 ) 
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X 


® 9 

= xq + x 0 h + h 


3 

y c f + o(h 5 ) 

L-i K K 


A 

X 


= x 0 + x„h + h 2 • Z“ 1[ f « + 0(h ' ) 

K = 0 


3 


X = X 0 + h • £c f +0(h 5 ) 

K=0 

c = c for k = 0, 1, 2 

K K 

c 3 = 0 

A 

C 4 = c 3 

Q! 3 = a 4 = 1 


( 71 ; 


(72; 


(73; 


Table 14 shows the reduced equations of condition for a fourth order 
formula RKN 4(5). 


21. The first four equations of Table 14 determine the weight factors and 

c as functions of the a t s. 
k K 


Since we require that the fifth evaluation should be taken over as first 
evaluation for the next step, the coefficients y 4x must be equal to the 

weight factors c (X = 0, 1, 2, 3). Setting y 31 = 0, the coefficients y 21 

and y 32 can then be determined from the two equations in the last line of 
Table 14. 


22. In Table 15 are listed the coefficients of a fourth-order Runge-Kutta- 

Ny strom formula RKN 4(5) and in Table 16, the error coefficients in x 
as well as in x for the formula of Table 15. There is only one error 
coefficient in x which is different from zero. The six error coefficients 
in x which are all different from zero are equal in magnitude or smaller 
than the error coefficient in x. Again, if the error coefficients differ 
by a constant factor only, those with the largest factor are listed in 
Table 16. 
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SECTION VI. SOME RUNGE-KUTTA-NYSTROM 
FORMULAS OF OTHER AUTHORS 

23. In Section VII we shall apply the Runge-Kutta-Nystrom formulas of 

Sections I through V to a numerical example. We shall compare their 
performance with that of our earlier Runge-Kutta formulas for first- 
order differential equations I 1] , 1 2] and also with the performance of 
some Runge-Kutta-Nystrom formulas by other authors. 

For the convenience of the reader we list in Tables 17 through 22 the 
Runge-Kutta-Nystrom coefficients and also the error coefficients for a 
fourth-, a fifth-, and a sixth-order Runge-Kutta-Nystrom formula, which 
we have used for comparison in Section VII. As before, if the error 
coefficients differ only by a constant factor, we list those with the largest 
factor. The fourth- and the fifth-order Runge-Kutta-Nystrom formulas 
were published by E. J. NYSTROM ([4], p. 24) and the sixth-order 
formula by J. ALBRECHT ([5], p. 103). 

A comparison of Table 18 with Table 16 shows that NYSTROM's error 

coefficients are two to nine times as large as ours. Therefore, 
NYSTROM's formula of Table 17 can be expected to have a larger lead- 
ing truncation error term than our formula RKN 4(5) . 

A comparison of Table 20 with Table 13 shows again that our formula 
RKN 5(6) has considerably smaller error coefficients than NYSTROM's 
formula RKN 5. For instance, our error coefficient T 5 is only about 
one -seventeenth of NYSTROM's coefficient T 5 . 

Again, the error coefficients of our Runge-Kutta-Nystrom formula 
RKN 6(7), as listed in Table 10, are smaller than those of ALBRECHT'S 
formula, listed in Table 22. For instance, our largest error coefficient 
in x, T 9 , is only about one-fifth of ALBRECHT'S coefficient T 9 . 

24. The formulas of NYSTROM and ALBRECHT of this section do not include 
an automatic stepsize control as our formulas of Sections I through V do. 
Therefore, we have to apply to the formulas of NYSTROM and ALBRECHT 
the standard stepsize control procedure that consists of recomputing two 
consecutive steps of stepsize h by one step of double stepsize 2h. It can 
easily be shown that the truncation error after one step of stepsize h 
is approximately 

E= -r • A x 2 , (74) 

2(2 - 1 ) 
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with n being the order of the formula under consideration and A x 2 the 
difference of our two results after two steps of stepsize h or one step of 
stepsize 2h, respectively. If the formula under consideration requires 
m evaluations per step, this stepsize control procedure would increase 
the number of evaluations per step to 2m - 1. So NYSTROM's formulas 
RKN 4 and RKN 5 and ALBRECHT'S formula RKN 6 would require five, 
seven, and nine evaluations per step, respectively. For NYSTROM's 
formulas, this is one evaluation more per step than our formulas 
RKN 4 (5) RKN 5(6) require, and for ALBRECHT’S formula RKN 6, 
this is two’ evaluations more than our formula RKN 617) requires, since 
for our formulas we can take over the last evaluation as first evaluation 
for the next step. Only for the very first integration step our formulas 
would require five, seven, and eight evaluations, respectively. Com- 
pared with NYSTROM's or ALBRECHT'S formulas, our formulas have 
the further advantage that they have considerably smaller leading trunca- 
tion error terms and therefore permit the use of a larger integration 
stepsize without loss of accuracy. The numerical example of Section VII 
will demonstrate this advantage. 


SECTION VII: APPLICATION TO A NUMERICAL PROBLEM 


25, In this section we apply the Runge-Kutta-Nystrb'm formulas of this report 
and some of the Runge-Kutta formulas of our earlier reports [ 1] , [2] to 
a numerical problem. For comparison, we also apply the Runge-Kutta- 
Nystrom formulas of Section VI to the same problem. 

In Table 23 the problem is stated, and the results of the numerical inte- 
gration are presented for the various formulas. 

All calculations were executed on an IBM-7094 computer in double preci- 
sion ( 16 decimal places) . The computer was equipped with an electronic 
clock to measure the running time for the various formulas. 


26. The stepsize control for our Runge-Kutta-Nystrom formulas RKN 4(5) , 
RKN 8(9) was set up in the following way. For a preset tolerance 
TOL (in Table 23 we used TOL = 0. 1 • 10“ 16 ), we computed for each 
step the products TOLX=TOL • |x 0 | , TOLY=TOL • ly 0 l which we con- 
sider as the tolerable errors in x and y for the step under considera- 
tion. Having computed the approximate truncation errors TE^, TE 
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according to Tables 3, 6, 9, 12, or 15, respectively, we then deter- 
mined the maximum of the ratios |TE f/TOLX, |TE l/TOLY and then 

x y 

required that for this maximum (max) the following inequalities hold: 
n + 1 

S max i 1 , (75) 



n being the order of our formula. 


If necessary, we halved or doubled the stepsize until (75) held. However, 
since the values TE^, TE^ are only approximations of the true trunca- 
tion error, it can happen that (75) never holds: If a certain stepsize h 

/ l \n + 1 

is too small max c I — J , the stepsize 2h might be too large: 


max > 1. In such a case, we accepted the smaller stepsize h as final. 
The stepsize control for our Runge-Kutta formulas RK 4(5) , . . . , 

RK 8(9) was set up quite similarly, but here we also tested the trunca- 
tion errors in x and y , since for these formulas we had to rewrite the 
differential equations of Table 23 as a system of four first-order dif- 
ferential equations in x, y, x, y . 


In the case of the Runge-Kutta-Nystrom formulas of Section VI we used 
equation (74) as approximation for the truncation error. 


27. Table 23 shows the result of the various formulas for t = 10. Since our 
problem has a solution in closed form, the total errors Ax, Ay, Ax, 

Ay are easily available and are listed in Table 23 for t = 10. In each of 
the Five groups of formulas, these errors do not differ much from one 
another. This means that all fourth-order formulas, etc. , are of about 
the same accuracy. 

However, the formulas of each group differ considerably from one 
another with respect to the number of steps required to cover the interval 
from t. = \j 7r/ 2 to t = 10 and with respect to the execution time on the 

computer. Our new Runge-Kutta-Nystrom formulas RKN 4(5) , . . , , 

RKN 8(9) require half or less than half the execution time of our earlier 
Runge-Kutta formulas RK 4(5) , . . . , RK 8(9) or of the formulas of 
Section VI. 


Computation Laboratory 

George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama 35812, September 9, 1971 
014-00-00-00-62 
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TABLE 1. COMPLETE EQUATIONS OF CONDITION FOR 
EIGHTH-ORDER FORMULA RKN 8( 9) 
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TABLE 1. COMPLETE EQUATIONS OF CONDITION FOR 
EIGHTH-ORDER FORMULA RKN 8(9) (Continued) 
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TABLE 1. COMPLETE EQUATIONS OF CONDITION FOR 
EIGHTH -ORDER FORMULA RKN 8(9) (Continued) 
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TABLE 1. COMPLETE EQUATIONS OF CONDITION FOR 
EIGHTH-ORDER FORMULA RKN 8(9) (Continued) 





TABLE 1. COMPLETE EQUATIONS OF CONDITION FOR 
EIGHTH-ORDER FORMULA RKN 8(9) (Continued) 





TABLE 1. COMPLETE EQUATIONS OF CONDITION FOR 
EIGHTH-ORDER FORMULA RKN 8(9) (Continued) 







TABLE 1. COMPLETE EQUATIONS OF CONDITION FOR 
EIGHTH-ORDER FORMULA RKN 8(9) (Concluded) 
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TABLE 3. COEFFICIENTS FOR AN EIGHTH-ORDER FORMULA RKN 8( 91 








TABLE 4. ERROR COEFFICIENTS FOR RKN 8(9) OF TABLE 3 





TABLE 5. REDUCED EQUATIONS OF CONDITION FOR 
SEVENTH -ORDER FORMULA RKN 7(8) 





TABLE 6. COEFFICIENTS FOR A SEVENTH-ORDER FORMULA RKN 7(8) 



RkSI 



TABLE 7. ERROR COEFFICIENTS FOR RKN 7(8) OF TABLE 6. 




TABLE 8. REDUCED EQUATIONS OF CONDITION FOR 
SIXTH-ORDER FORMULA RKN 6(7) 
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12 

= C 2 «| + c 3 «| + C 4 cy| + c 5 + c 6 
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j_ 

20 

= C 2 a\ + c 3 o-l + c 4 QfJ + c 5 a\ + c 6 
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J_ 

42 
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= c 2 P 23 + c 3 P 33 + c 4 P 43 + c 5 P 53 + c 6 • 20 
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" 120 

h g 
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TABLE 9. COEFFICIENTS FOR A SIXTH-ORDER FORMULA RKN 6 ( 7 ) 
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TABLE 10. ERROR COEFFICIENTS FOR RKN 6(7) OF TABLE 9 
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TABLE 13. ERROR COEFFICIENTS FOR RKN 5(6) OF TABLE 12 



TABLE 14. REDUCED EQUATIONS OF CONDITION FOR 
FOURTH-ORDER FORMULA RKN 4( 5) 
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TABLE 15. COEFFICIENTS FOR A FOURTH-ORDER FORMULA RKN 4(5) 



TABLE 16. ERROR COEFFICIENTS FOR RKN 4(5) OF TABLE 15 
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TABLE 17. COEFFICIENTS FOR NYSTROM'S 
FOURTH-ORDER FORMULA RKN 4 



TABLE 18. ERROR COEFFICIENTS FOR NYSTROM'S 
FOURTH-ORDER FORMULA RKN 4 OF TABLE 17 

















TABLE 19. COEFFICIENTS FOR NYSTROM'S 
FIFTH-ORDER FORMULA RKN 5 



TABLE 20. ERROR COEFFICIENTS FOR NYSTROM'S 
FIFTH -ORDER FORMULA RKN 5 OF TABLE 19 

T 2 ^ \ ( <-T 0 i\ + c 2 a\) - 7^ ~ " 0- °°0 277 778 
T 5 r 7 c 2 P 22 - ~ as - o. 000 555 556 

r 3 | ( <7 cv? + c 2 cv 2 + c 3 of§) - 000 158 889 

• 1 . . 1 

[ G 7 U 2 or 2 P 22 + c 3 <* 3 P 32 ) - — as 0. 000 555 556 
~ 144 

- 7 (6 2 p 23 + C3 P 33> - 7^7 w 0. 000 277 778 







TABLE 21. COEFFICIENTS FOR ALBRECHT'S 
SIXTH-ORDER FORMULA RKN 6 



59 




























TABLE 22. ERROR COEFFICIENTS FOR ALBRECHT'S 
SIXTH-ORDER FORMULA RKN 6 OF TABLE 21 

! T 3 ( Cj + t’2 a 2 + c 3 Off J - 7~7T « - o. 000 04G 503 _ 

c 4 oob 

1 o o 1 

T 4 — (c 2 a\ P 21 + c 3 a 3 P 31 ) - — » - 0.000 074 405 

•*■6 7 ( c 2 a 2 f*22 + c 3 a 3 ^32^ - ~ qq 8 * - 0. 000 037 202 

T a 7 < °2 p 23 + c 3 P 33 i - - * gQ » - 0. 000 074 405 

T 10 = C 3 y 32 P 21 — as 0. 000 032 873 

5 040 

'i’a = ( Cl O!® + C 2 Of® + c 3 a® + c 4 a\) - « 0. 000 023 251 

Tg 7 (c 2 a 2 P 21 + c 3 aj P 31 + c 4 aj P 41 ) - ~ « 0. 000 074 405 

f 7 - ~ ( c 2 a\ p 22 + c 3 a l P 32 + c 4 arf P 42 ) - ^ » 0. 000 018 601 

1 

T 9 7 (c 2 P 21 + c 3 Pf< + c 4 P| 4 ) - YYY * °- 000 066 138 

3 13 ■= 7 (c 2 Of 2 P 23 + C 3 Of 3 P 33 + c 4 a 4 P 43 ) - « 0. 000 074 405 

P 14 c 3 c: 3 y 32 P 21 + c 4 a 4 (y 42 P 21 + y 43 P 31 ) - ~ - 0. 000 033 069 

p 16 - \ ( P 24 + c 3 P 34 + c 4 P 44 I ~ g^Q * 0. 000 046 503 

P 18 C 3 • y 32 a 2 P 21 + c 4 (y 42 a 2 P 24 + y 43 a 3 P 3 i) - — " - g Q « 0 . 000 074 405 

j 19 “ TJ ^ 3 ^ 32 P 22 + ^4 (Y42 P 22 + y 43 P 32 ) ] “ 5~ 040 ~ ^18 8 ^1 
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TABLE 23. APPLICATION OF THE VARIOUS FORMULAS TO A NUMERICAL PROBLEM 











